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Abstract 

The "developmental hourglass" describes a pattern of increasing morphological 
divergence towards earlier and later embryonic development, separated by a 
period of significant conservation across distant species (the "phylotypic 
stage"). Recent studies have found evidence in support of the hourglass effect 
at the genomic level. For instance, the phylotypic stage expresses the oldest 
and most conserved transcriptomes. However, the regulatory mechanism that 
causes the hourglass pattern remains an open question. Here, we use an 
evolutionary model of regulatory gene interactions during development to 
identify the conditions under which the hourglass effect can emerge in a 
general setting. The model focuses on the hierarchical gene regulatory network 
that controls the developmental process, and on the evolution of a population 
under random perturbations in the structure of that network. The model 
predicts, under fairly general assumptions, the emergence of an hourglass 
pattern in the structure of a temporal representation of the underlying gene 
regulatory network. The evolutionary age of the corresponding genes also 
follows an hourglass pattern, with the oldest genes concentrated at the 
hourglass waist. The key behind the hourglass effect is that developmental 
regulators should have an increasingly specific function as development 
progresses. Analysis of developmental gene expression profiles from 
Drosophila melanogaster and Arabidopsis thaliana provide consistent results 
with our theoretical predictions. 
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Introduction 

The evolutionary mechanism of conservation during embryogen- 
esis, and its connection to the gene regulatory networks that control 
development, are fundamental questions in systems biology 13 . Sev- 
eral models have been presented in the context of morphological, 
molecular, and genetic developmental patterns. The most widely 
discussed model is the "developmental hourglass", which places 
the strongest conservation across species in the "phylotypic stage". 
The first observations supporting the hourglass model go back to 
von Baer when he noticed that there exists a mid-developmental 
stage in which embryos of different animals look similar 4 . On the 
other hand, the "developmental funnel" model of conservation pre- 
dicts increasing diversification as development progresses 5 6 . 

Recently, the hourglass model has come under new light. Multiple 
studies have observed the hourglass pattern across diverse biologi- 
cal processes, including transcriptome divergence 7-11 , transcrip- 
tome age 71213 , molecular interaction 14 , and evolutionary selective 
constraints 101415 . Despite these observations the genomic basis and 
even the existence of the developmental hourglass effect have been 
the subject of an intense debate 1 ' 61316-21 . More importantly, the under- 
lying mechanism that can shape the developmental process in the 
hourglass or funnel forms is still unknown. 

We aim to understand the conditions under which the hourglass 
effect can emerge in a general setting, based on an abstract model 
for the evolution of embryonic development. The model focuses 
on a hierarchical network that represents the temporal "execution" 
of the underlying Gene Regulatory Network (GRN) during devel- 
opment. Each layer of the network corresponds to a developmen- 
tal stage. The nodes at each layer represent regulatory genes (i.e., 
genes encoding transcription factors or signaling molecules) that 
undergo significant activity change at that corresponding stage. The 
edges from genes at one layer to genes at the next layer represent 
regulatory interactions that cause those activity changes. We refer 
to this hierarchical network as Developmental Gene Execution 
Network (DGEN) to distinguish it from the corresponding GRN. 
A DGEN is subject to evolutionary perturbations (e.g., gene dele- 
tions, rewiring, duplication) that may be lethal, or that may impede 
development, for the corresponding organism. 



The model predicts that the evolutionary process shapes the DGENs 
of a population in the form of an hourglass, under fairly general 
assumptions. Specifically, the number of genes at each developmen- 
tal stage follows an hour-glass pattern, with the smallest number at 
the "waist" of the hourglass. The main condition for the appearance of 
the hourglass pattern is that the DGEN should gradually get sparser as 
development progresses, with general-purpose regulatory genes at 
the earlier developmental stages and highly specialized regulatory 
genes at the later stages. Another model prediction is that the evolu- 
tionary age of DGEN genes also follows an hourglass pattern, with 
the oldest genes concentrated at the waist. 

We have examined the aforementioned predictions using transcrip- 
tome data from the development of Drosophila melanogaster and 
Arabidopsis thaliana. This data is insufficient to reconstruct the 
complete DGEN of these species but it allows to estimate the num- 
ber of genes at each developmental stage, given an activity variation 
threshold. Under a wide range of this threshold, the inferred DGEN 
shape follows an hourglass pattern, the waist of that hourglass 
roughly coincides with the previously reported phylotypic stage for 
these species, and the age of the corresponding genes follows the 
predicted hourglass pattern. 

Developmental gene execution networks 

As a first-order approximation, a regulatory gene can be modeled 
in one of several discrete functional states 12 . In the simplest case, a 
regulatory gene can act as a binary switch ("on" and "off) but in 
general a gene may have more than two functional states. The tran- 
sition of a regulatory gene X from one functional state to another is 
often (but not always) caused by one or more upstream regulators 
of X that go through a functional state transition before X. We use 
the term transitioning gene to refer to a regulatory gene that goes 
through a functional state transition at a given developmental time 
anywhere at the developing embryo. 

A DGEN is a directed and acyclic network; see Figure la for an 
abstract example. The vertical direction refers to developmental 
time, from the zygote at the top to the developed organism at the 
bottom. In the horizontal direction we can represent different spa- 
tial domains, even though this is not necessary and it is not done in 




DF 



Figure 1. An abstract DGEN. (A) The circles denote state-transitioning genes, edges represent directed regulatory interactions, and colored 
boxes refer to spatial domains that form during development. If regulatory genes become increasingly function-specific as development 
progresses, the network gradually becomes sparser in that direction. (B) Evolutionary perturbations on a DGEN's structure: Gene A is deleted 
(DL), while gene B is rewired (RW) losing an outgoing edge. This RW event causes the regulatory failure (RF) of gene C, which then causes 
a cascade of five more RF events. This cascade causes developmental failure (DF). Note that the removal of some upstream regulators does 
not always cause an RF event (e.g., genes regulated by A). 
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our model. For instance, the zygote at the top of the DGEN would 
be a single domain, while the organism at the bottom of the DGEN 
would have the largest number of spatial domains. Development is 
often approximated (conceptually and experimentally) as a succes- 
sion of discrete developmental stages. The duration of a develop- 
mental stage can be thought of as the typical time that is required 
for a gene's functional state transition, and it does not need to be the 
same for all stages. Each layer of a DGEN refers to a developmen- 
tal stage, and it includes only the transitioning genes during that 
stage anywhere in the embryo. The same gene can appear in more 
than one stage if it goes through several functional state transitions 
during development. Additionally, a DGEN edge from a gene X at 
stage / to a gene Y at stage / + 1 implies that the functional transition 
of X caused the functional transition of Y at the next stage. If gene Y 
has more than one incoming edge, its functional state transition was 
caused by the coupled effect of more than one transitioning genes 
at the earlier stage. Any upstream regulators of Y that remained at 
the same functional state at stage / are not included in that stage of 
the DGEN. 

The sequence of developmental stages is denoted by 1=1 ... L. The 
set of transitioning genes at stage / is G(/). A gene g at stage 1<L 
regulates a set of downstream genes at stage / + 1 denoted by D(g) 
(outgoing edges from g). Similarly, a gene g at stage />1 is regu- 
lated by a set of upstream genes U(g) at stage / - 1 (incoming edges 
to g). The functional transitions at the first stage are assumed to be 
caused by regulatory maternal genes that initiate the developmental 
process. 

Model description 

The model captures certain aspects of both the developmental pro- 
cess, in the form of a given DGEN for each embryo, and of the 
evolutionary process, as random perturbations in the structure of 
individual DGENs in a population. The model does not need to cap- 
ture the actual functional state transitions or the regulatory input 
function of each gene. It does capture however the dynamic and 
stochastic effect of structural network perturbations (gene deletion, 
rewiring and duplication) on the success of the developmental pro- 
cess, as explained in the following. 

Similar to the Wright-Fisher model, we consider a population of 
N individuals, each represented by a DGEN. In each generation, 
individuals reproduce asexually, inheriting the DGEN of their par- 
ent. Various evolutionary events can cause structural changes in the 
DGEN of an individual that may result in "developmental failure". 
Such individuals (and their DGENs) are replaced with develop- 
mentally successful individuals so that the population size remains 
constant. 

The model is meant to be as general as possible and so the regula- 
tory interactions between genes of successive stages are determined 
probabilistically, as follows. Each stage / is assigned a regulatory 
specificity, or simply specificity s(l) with 0 < s(l) < 1. A gene g at 
stage / acts as upstream regulator for a gene g' at stage / + 1 with 
probability s'{l) = 1 - s(l). So, the specificity of a developmental 
stage determines how likely it is for regulatory genes of that stage 
to cause a state transition of the genes at the next stage; a higher 
specificity decreases that likelihood. 



Our major assumption is that the regulatory specificity increases 
substantially as development progresses. In other words, the DGEN 
becomes gradually sparser along the developmental time axis, start- 
ing with s(l)»0 and ending with s(L)«l. This assumption is plausible 
for the following reasons. First, as development progresses the embryo 
grows in size forming distinct spatial domains. So, extracellular 
gene regulation becomes more difficult, especially across different 
domains. Additionally, as development progresses the transitioning 
genes become more organ- or tissue-specific, implying that their 
downstream interactions become sparser. Unfortunately, an empiri- 
cal investigation of the increasing specificity assumption requires 
knowledge of the complete DGEN for a given species; this is cur- 
rently not feasible for even the most well-studied model organisms. 

The DGEN structural changes we consider are gene deletions, gene 
duplications, and gene rewiring: 

Deletions (DL): This event removes a gene from the DGEN, includ- 
ing its incoming and outgoing edges. There are many genetic mech- 
anisms that may cause such events. A DL event deletes each gene of 
an individual and at each generation with probability P DL . 

Duplication (DP): This event creates an identical copy of a gene g 
with the same downstream and upstream regulators and at the same 
developmental stage as g. The two genes may have different fates 
if one of them is subject later to deletion or rewiring. Otherwise, 
the two genes are considered identical. A DP event duplicates each 
gene of an individual and at each generation with probability P Dp . 

Rewiring (RW): This event changes the upstream and/or down- 
stream regulators of a gene. Changes in the upstream versus down- 
stream regulators may have different biological basis. The former 
occur, for instance, as a result of mutations in the transcription factor 
binding sites in a gene's promoter or mutations in distal regulatory 
elements such as enhancers, while the latter may be mostly caused 
by coding sequence mutations. The details of the rewiring process 
do not affect the results qualitatively as long as the average density 
of edges in each stage remains consistent with the specificity of that 
stage. The specific rewiring mechanisms we use are presented next. 

Suppose that a RW event affects gene g at stage /. The upstream 
regulators of g are recomputed based on the specificity of the previ- 
ous stage, i.e., by choosing each distinct gene of stage / - 1 with 
probability s f (l - 1). For the downstream regulators of g, we ran- 
domly remove N_ existing outgoing edges of g, and then add A/+ 
outgoing edges to randomly chosen genes of stage / + 1 that g is not 
already connected to. Both N_ and N+ follow a Binomial distribution 
with \D(g)\ trials and success probability s f (l). This captures that 
the downstream regulators of g are derived by incremental changes 
in D(g), instead of giving g a completely new network configura- 
tion (thereby, new regulatory function). The higher the regulatory 
specificity of a stage, the less likely these incremental changes are. 
An RW event rewires each gene of an individual and at each genera- 
tion with probability P RW . 

A gene deletion or rewiring event at stage / can remove an upstream 
regulator from genes at stage / + 1. A loss of incoming edges may 
trigger the regulatory failure of a gene, as described next. 
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Regulatory failures (RF): A gene g may not be able to change 
functional state if some of its upstream regulators U(g) are lost 
due to DL or RW events (see Figure lb). Even though regulatory 
networks are often robust to structural perturbations, even a partial 
gene loss in U(g) may disable g causing a regulatory failure. It is 
plausible that the probability of a regulatory failure increases with 
the fraction of lost upstream regulators. So, if U'(g) is the new set 
of upstream regulators and \U(g)\ > \U'(g)\ > 0, gene g is removed 
with probability: 



PRF(r) = \- 



A-r 



0<r = l- 



\U'(g)\ 
P(g)\ 



<1 



while if \U'(g)\ = 0 we set P RF (1)=1. z is the RF parameter and it 
depends on the robustness of regulatory interactions to gene loss 
(see Figure 2). 

When a DL or RW event causes one or more RF events, the latter can 
trigger additional RF events in subsequent developmental stages, 
leading to cascades of regulatory failures. Such RF cascades may 
cause developmental failure, meaning that the developed embryo is 
unable to survive or reproduce. 



Methods 

Hourglass score H. Suppose that w(l) denotes the number of tran- 
sitioning genes in stage / and let b be the stage with the minimum 
number of such genes. We construct the sequences X={w(l)}, I = 
1, ... b) and Y={w(l)}, I = b, ... L}. t x and T y denote the normal- 
ized univariate Mann- Kendall statistic for monotonic trend in each 
sequence, respectively (t is -1 for decreasing, +1 for increasing 
and almost 0 for random sequences). The H score is defined as 
H = (t y - Ty)/2. See Figure 3 for an illustration, and for the defini- 
tion of a more robust version of H. 

Statistical analyses. All biological data processing and analyses 
were performed using custom scripts written in Java (JDK vl.6) 
[Dataset. 14]. 

Drosophila data and treatment. Developmental gene expression 
profiles for D. melanogaster are obtained from two sources. First, 
we obtained microarray data from Kalinka et aV for 3,610 genes. 
The expression level of each gene is calculated as the median of 
probes mapping to that gene. Each stage represents a 2-hr interval 
during the first 20 hours of embryogenesis (10 stages). The sec- 
ond source is RNA-Seq data from Graveley et alP. Raw data are 



Developmental failure (DF): The last stage of a DGEN represents 
the fully developed embryo. If that stage includes T transitioning 
genes at a successfully developed embryo, the simplest assumption 
is that an individual with less than T genes at stage-L has failed to 
develop properly. Such DGENs are removed from the population 
and they are replaced with randomly chosen but successfully devel- 
oped DGENs. We have also experimented with two variations of 
the DF event: first, the individual is removed if its last stage has 
less than T - /genes, where yis small relative to T, and second, the 
probability of a DF event increases as the number of genes at stage-L 
decreases below T. The qualitative results, as described next, do not 
change with these two model variations. 




0 0.2 0.4 0.6 0.8 



r 

Figure 2. Probability of Regulatory Failure (RF) for three values 
of the paramater z. r is the fraction of upstream regulating edges 
that are lost due to a DL or RW event. 



Stage width 



w(l) 




Figure 3. Illustration of the H score calculation. Let w{l) be the 
width of stage /. Let w b be the minimum width across all stages, and 
suppose that this minimum occurs at stage / = b; this is the waist 
of the network (ties are broken so that the waist is closer to lU2\). 
Consider the sequence X = {w{l)}, I = 1 , ... b) and the sequence Y 
= {w(l)}, I = b, ... /_}. We denote the normalized univariate Mann- 
Kendall statistic for monotonic trend on the sequences Xand Vas 
r x and T r respectively. The Mann-Kendall statistic varies between 
-1 (decreasing) and 1 (increasing), and it is approximately zero 
for a random sequence. We define H = (t y - r x )/2; H is referred 
to as the hourglass score. H = 1 if the DGEN is structured as an 
hourglass, with a decreasing sequence of b stages followed by an 
increasing sequence of L - b stages. In the computational modeling 
results, we do not consider the width of the first stage because it 
can never decrease in Models-1 ,2,3. We also define a variation of 
the hourglass score in which we do not take into account adjacent 
stages in calculating the two Mann-Kendall statistics. That statistic 
is denoted by Hand is referred to as the robust hourglass score. 
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processed to RPKM values. Each stage represents a 2-hr interval 
during the first 24 hours (12 stages). Genes with zero RPKM value 
in all stages are discarded, resulting in 14,110 genes. 

Arabidopsis data and treatment. Genome-wide expression profiles 
of a complete developmental series from the zygote to the mature 
embryo in A. thaliana were obtained from Xiang et al. 24 . This com- 
prised of microarray expression levels for 25,207 genes across seven 
developmental stages. Signal background correction and normaliza- 
tion of raw expression values was carried out by Xiang et al. 24 . 



The transcriptome age index (TAI) values for Arabidopsis genes 
were obtained from 11 . 

Age index for each stage-pair. Suppose that we identify transi- 
tioning genes based on the normalized expression levels, and that 
n(l) genes are assigned to stage-pair (/ - 1, /). Denote by p. the 
phylogenetic rank (TAI value) of gene i. The age index assigned to 
that stage-pair is given by TAl(l)=(H'^p.e'.^)/ll'^p.e'. r The same 
method is used when transitioning genes are identified based on 
absolute expression levels. 



Transitioning gene identification. Suppose that the reported expres- 
sion value of gene i at stage / is e. f We analyze both these "abso- 
lute" expression values as well as the normalized expression values, 
given by e' . =e.,/H e.,. 

The identification of transitioning genes follows the same method 
for both absolute and normalized expression levels. In the case of 
normalized expressions, we calculate 8 n = e' u - e'. l _ l for each gene 
and at each stage 1=2 ... L. Gene i is considered "transitioning" 
at the stage-pair (/ - 1, /) if |<5.J > c, where c is a given transition 
threshold. This condition is more robust to noise than a ratio-based 
rule (e'./e'.^) for the identification of transitioning genes. Note 
that a gene may be transitioning at more than one stage-pair, but it 
may also not be transitioning at any stage-pair. 

Transcriptome age index (TAI). We collected the groups of orthologs 
for each gene in Drosophila using two databases, OrthoDB 25 and 
OrthoMCL 26 . The Eumetazoa data were taken from OrthoDB, while 
Fungi and Plants species were retrieved from OrthoMCL, and the two 
datasets were merged. Using these orthologs we then assigned an age 
index to each gene based on its absence and presence in a phylogenetic 
tree of 24 well-diverged species 1112 [Dataset. 7] (see Figure 4). 



Results 

Simulation 

We simulate the presented model to examine the properties of the 
surviving DGENs as evolutionary time progresses. The initial pop- 
ulation consists of N identical DGENs with T genes at each stage. 
The edges between genes are constructed probabilistically based on 
the specificity of each stage, as described previously. Simulating 
the complete model would not show the significance of individual 
mechanisms such as the increasing specificity assumption. For this 
reason we construct a sequence of four models with increasing 
complexity, presenting results separately for each of them: 

Model- 1: Constant specificity. Each stage has the same specific- 
ity, s(l)=0.5 for / = 1 ... L - 1. Further, this model does not include 
gene deletion and duplication. Gene rewiring can cause RF and DF 
events even if there are no DL or DP events. 

Model-2: Increasing specificity. The difference from Model- 1 is 
that the specificity is gradually increasing across developmental 
stages. Unless noted otherwise, the specificity is linear, s(l) = l/L 
for /= 1 ... (L - 1); a nonlinear specificity function is also consid- 
ered, which we describe later. 




D melanogaster 
D pseudaabscura 
A gambiae 
N vitripennis 
A mellifera 
D pulex 
I scapularis 

M rnusculus 
H sapiens 
G gallus 
Xtropicalis 
T nigrovi ridis 
S purpuratus 
C elegans 
C briggsae 

vectensis .„ 
magnipapillata 

N crassa 
A niger 
Y lipolytica 
S cerevisiae 
E cuniculi 
A thaliana 
0 sativa 



Figure 4. The phylotypic tree that we use to calculate the age index of Drosophila's genes. Each gene is assigned to one of the following 
six ages: Level- 1 : Common ancestor to Fungi, Plants and Eumetazoa. Level-2: Common ancestor to Fungi and Eumetazoa. Level-3: Common 
ancestor to all Eumetazoa. Level-4: Common ancestor to all Bilateria. Level-5: Common ancestor to all Arthropoda. Level-6: Common ancestor 
to all Dipteria. 
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Model-3: With gene duplications. Model- 3 adds DP events in 
Model-2. The duplication probability P Dp is set so that the average 
size of a DGEN, across the entire population, stays within a given 
range (70%-80% ofLxT genes). 

Model-4: With gene deletions. Model-4 adds DL events in Model-3 
(complete model). The deletion probability P DL is set so that the 
average size of a DGEN, across the entire population, stays within 
the same range as in Model-3. 

In Model- 1 and Model-2, genes can be only removed (due to RW 
events, potentially followed by RF cascades) and so the average 
DGEN size decreases as evolutionary time progresses, which is 
unrealistic. Model-3 and Model-4 are more realistic because they 
can maintain a roughly constant DGEN size in the long-term. 
However, as will be shown next, all aspects of the developmental 
hourglass effect can already be seen with Model-2 (but not with 
Model- 1). This highlights the increasing specificity assumption as 
the key property behind the developmental hourglass effect. 

Hourglass shape. A first observation is that as evolutionary time pro- 
gresses, DGENs acquire an "hourglass-like" shape in Models-2,3,4. 
This means that the width of each stage first decreases until a 



certain stage (referred to as the waist of the hourglass) and then grad- 
ually increases. The hourglass may not be symmetric with respect 
to the waist. To quantify this observation, we define an "hourglass 
score" H (see Methods and Figure 3) that is equal to 1 if the sequence 
of L stage widths consists of two segments: a decreasing sub- 
sequence of k > 2 stages followed by an increasing sub- sequence of 
L - k > 2 stages. Figure 6a shows the hourglass score for the popula- 
tion of DGENs in Model-2. Similar graphs for the three other mod- 
els are shown in Figure 5a, Figure 7a, and Figure 8a. The H score 
quickly increases in the three models that have increasing specific- 
ity, and it fluctuates close to 1 afterwards. What is the reason behind 
the hourglass shape of DGENs? When a gene g is rewired at stage /, 
it may trigger RF events in stage / + 1 depending on the number of 
its lost outgoing edges. In the first few stages, where specificity is 
low, it is unlikely that a gene would lose a large fraction of its (typi- 
cally many) incoming edges. In the last few stages, where specificity 
is high, edges are unlikely to get rewired in the first place. In the mid- 
stages however, where the specificity is close to 50%, there is higher 
variability in the number of outgoing edges lost or gained due to 
RW events. The loss of several outgoing edges due to an RW event 
at stage / can trigger several RF events and gene removals in the 
subsequent stage. Thus, the probability of RF events in mid-stages is 
higher than in early/late stages, making the removal of genes more 
likely in the former. 




600000 800000 





Figure 5. Computational results with Model-1. Parameters: 10 runs with different initial populations, N=^0 individuals, L=^0 stages, 
specificity function s(/)=0.5 for all stages, r=100 genes at each stage initially, RF parameter z=4, 1,000,000 generations, probability of RW 
event P RW =10~ 4 . The red line is the median and the green boxes are the 10th, 25th, 75th, and 90th percentiles, across all individuals and 
all runs, (a) The hourglass score H across evolutionary time, (b) Lethality probability at each stage, (c) Age of existing genes at the last 
generation. 




Generation Stage Stage 

Figure 6. Computational results with Model-2. Parameters: 10 runs with different initial populations, A/=1000 individuals, L=^0 stages, 
specificity function s( l)=l/L (/=1 ... L - 1), r=100 genes at each stage initially, RF parameter z=4, 500,000 generations, probability of RW event 
P flM = 1fj- 4 . The red line is the median and the green boxes are the 10th, 25th, 75th, and 90th percentiles, across all individuals and all runs. 
(A) The hourglass score /-/across evolutionary time. (B) Lethality probability at each stage. (C) Age of existing genes at the last generation. 
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Figure 7. Computational results with Model-3. Parameters: 10 runs with different initial populations, N=^0 individuals, L=^0 stages, 
specificity function s(l)=l/L (/=1 ... L - 1), r=100 genes at each stage initially, RF parameter z=4, 1,000,000 generations, probability of RW 
event P RW = 10 -4 . The probability of gene duplication P DP is adjusted dynamically so that the average DGEN size stays between 700 and 800 
genes. The red line is the median and the green boxes are the 10th, 25th, 75th, and 90th percentiles, across all individuals and all runs, (a) 
The hourglass score /-/across evolutionary time, (b) Lethality probability at each stage, (c) Age of existing genes at the last generation. 
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Generation 





Figure 8. Computational results with Model-4. Parameters: 10 runs with different initial populations, A/=10 individuals, /_=10 stages, 
specificity function s( /)=///_ (/=1 ... L - 1), r=100 genes at each stage initially, RF parameter z=4, 1,000,000 generations, probability of RW 
event P RW = 10~ 4 . The probability of gene duplication P DP is adjusted dynamically so that the average DGEN size stays between 700 and 800 
genes. The probability of gene deletion (DL) is P DL = 10~ 6 . The red line is the median and the green boxes are the 10th, 25th, 75th, and 90th 
percentiles, across all individuals and all runs, (a) The hourglass score /-/across evolutionary time, (b) Lethality probability at each stage, (c) 
Age of existing genes at the last generation. 



The constant specificity of Model- 1 does not result in an hourglass 
pattern [Dataset. 1] (see Figure 5a) for the following reason. RW 
events at stage / can cause RF events at the next stage with the same 
probability, independent of /. However, after the occurrence of an 
RF event, the size of the potential cascade increases as / decreases 
simply because there are more subsequent stages to affect. This 
gives DGENs a "funnel-like" shape with a gradually increasing 
number of genes after stage- 1; //fluctuates around 0.5, as expected 
for an increasing sequence. 

Stage lethality. Another aspect of the developmental hourglass is 
in terms of the significance of each stage for the survival of the 
embryo. We define lethality of stage I as the probability that a RW 
or DL event at stage / starts a RF cascade that eventually leads to a 
DF event. We estimate this probability at generation i as the fraction 
of RW and DL events, during the past i generations, that occurred at 
stage / and led to a DF. 

In Model- 1 , there is no clear trend for the stage lethality probabil- 
ity (see Figure 5b); with the exception of the last stage (in which 
RW events cannot result in gene loss), the lethality probability is 
roughly the same at all stages. In the three models with increasing 
specificity, however, we observe a clear pattern: the lethality gradu- 
ally increases until the waist of the hourglass, and then it decreases 
[Dataset. 2, Dataset. 3 and Dataset. 4] (see Figure 6b, Figure 7b, 



and Figure 8b). The reason, as explained earlier, is that the prob- 
ability of RF events in mid-stages is higher that in early/late stages. 
Additionally, after the formation of the hourglass shape the mid- 
stages have relatively few genes and so an RF event in those genes 
is more likely to initiate a lethal RF cascade. 

Age of genes. A third aspect of the developmental hourglass effect 
is related to the evolutionary age of genes. The age of a gene g at 
generation i is defined as A(g) = i - t Q (g), where t Q (g) is the gen- 
eration at which g was most recently rewired (and 0, if it was not 
rewired earlier). The rationale behind this definition is that a rewir- 
ing event may give that gene a new function, at least in terms of its 
upstream and downstream regulators. 

In the case of Model-2, Figure 6c shows the median age of the 
genes at each stage, considering the population of all genes across 
all individuals at a given generation. See Figure 5c, Figure 7c, and 
Figure 8c for the same results with the three other models. The evo- 
lutionary age at stage / follows the same pattern as the lethality 
probability: it gradually increases until we reach the waist of the 
hourglass, and then it gradually decreases. Genes at intermediate 
stages tend to be older because, as explained earlier, they are fewer 
and their rewiring is more likely to be lethal. When one of those 
genes g is rewired or deleted from a DGEN, the corresponding 
individual is often replaced (DF event) by another individual that 
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has the same gene g. So, the genes at the waist of a DGEN tend to 
be more conserved than genes at earlier or later stages. 

Location of waist. What controls the location of the hourglass waist 
in the developmental process? 

The location of the waist is mostly affected by two parameters of the 
model: the shape of the specificity function and the RF parameter z. 
To examine the effect of the former we use the nonlinear function 
shown in Figure 9. /is the stage at which the specificity is 50%, and 
so that stage has the maximum variance in the number of outgo- 
ing regulatory edges. RW events at this stage can cause the largest 
extent of rewiring and so, the highest likelihood of RFs in genes of 
the next stage. Figure 10a shows that the location of the hourglass 
waist is "pushed" towards stage y, even though it does not coincide 
always with that stage. Parameter z controls the shape of the RF 
probability: increasing z makes RF events more likely, also increas- 
ing the likelihood of lethal RF cascades. Figure 10b shows that as 
we increase z the hourglass waist moves towards later developmen- 
tal stages [Dataset. 5]. 

Gene prevalence. We introduce a "gene prevalence" metric for 
gene g at time t as the fraction of individuals that include g at evo- 
lutionary time t [Dataset. 6]. Figure 11a shows the gene prevalence 
metric across developmental stages after 500,000 generations, 
while Figure lib shows the relation between gene prevalence and 
gene age. The genes at the waist of the developmental hourglass are 
not only the oldest but also the most prevalent across the popula- 
tion. The implication of this simulation result is that we can expect 
that those genes that are transitioning near the waist of the develop- 
mental hourglass will be the most conserved genes in a population. 



Data analysis 

We have examined the predictions of the previous model using 
transcriptome data for Drosophila melanogaster and Arabidopsis 
thaliana [Dataset. 14]. We summarize results from both species 
here; the corresponding figures for which are Figure 12 and Figure 13 
[Dataset. 8 and Dataset. 9]. Vox Drosophila, we analyze Microarray 9 
and RNA-Seq 23 temporal expression profiles during the first 20 hours 
of development, taken at 10 stages of 2-hr intervals. We examine 
whether a) the number of transitioning genes follows an hourglass 
pattern, b) the waist of that hourglass coincides with the Drosophila 
phylotypic stage, and c) the evolutionary age of the transitioning 
genes follows a similar hourglass pattern. The two datasets are 
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Figure 9. A nonlinear specificity function, s(/)=0.9-^^- , for 

three values of the parameter y. This function allows us to control 
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Figure 10. We examine the effect of the two model parameters that affect the location of the DGEN hourglass waist. The first is the 
specificity function. To examine its effect, we use a sigmoid-like mathematical function that controls the stage y at which the specificity is 
50% (see Figure 9). This is the stage with the maximum variance in the number of outgoing regulatory edges. RW events at this stage can 
cause the largest extent of rewiring and so, the highest likelihood of RFs in genes of the next stage. Graph (a) shows that the location of the 
hourglass waist is "pushed" towards stage 7, even though it is not always exactly at that stage. The second way to affect the location of the 
hourglass waist is the parameter zthat controls the shape of the RF probability Increasing z makes RF events more likely, also increasing 
the likelihood of lethal RF cascades. Graph (b) shows that as we increase zthe hourglass waist moves towards later developmental stages. 
These results are obtained using Model-2. Parameters: 10 runs with different initial populations, A/=1000 individuals, Z_=10 stages, specificity 
function s(l)=l/L (/=1 ... L- 1), r=100 genes at each stage initially, RF parameter z=4, 500,000 generations, probability of RW event P RW = 
10" 4 . The graphs show the median (red lines) and the 10th, 25th, 75th, and 90th percentiles (green boxes) of the location of the waist. 
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described in more detail in the Methods section. With such limited 
data, we cannot infer the regulatory edges between transitioning 
genes and we cannot reconstruct the underlying DGEN. However, 
we can identify the transitioning genes at each developmental stage 
given a "transition threshold" c (see Methods). Even though the 



correct value of this threshold is not known, the following results 
are robust in a wide interval of c, which includes most of the expres- 
sion variation range across successive developmental stages (see 
Figure S3 for the CDFs of expression level variations across suc- 
cessive stages [Dataset. 12]). 





0.4 0.5 0.6 
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Figure 11. The prevalence of a gene g in a population of N individuals is the fraction of individuals in which gene g appears. These 
results are obtained using Model-2. Parameters: 10 runs with different initial populations, A/=1000 individuals, /_=10 stages, specificity function 
s(l)=l/L (/=1 ... Z_-1), r=100 genes at each stage initially, RF parameter z=4, 500,000 generations, probability of RW event P RW = 10 -4 . The 
graphs show the median (red lines) and the 10th, 25th, 75th, and 90th percentiles (green boxes) for: (a) prevalence of genes in each stage after 
500,000 generations, and (b) gene age as a function of gene prevalence. As expected, older genes tend to be more prevalent in the population. 
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Figure 12. Drosophila results using normalized expression levels. Graphs (A) and (B) show the hourglass score (normal and robust) 
as a function of the transition threshold c for the two datasets. Graphs (C) and (D) show the location of the hourglass waist (stage-pair) as 
a function of the transition threshold c for the two datasets. Graphs (E) and (F) show the Transcriptome Age Index of transitioning genes 
for three different values of c (chosen so that the number of genes with known age index assigned to each stage is at least 25) for the two 
datasets. Graph (G) shows the transitioning genes for the Microarray dataset with c=0.0005. The transitioning genes constitute 11% of all 
genes in that dataset. 53% of those genes transition in a single stage-pair. Of the remaining, 64% transition only in consecutive stage-pairs. 
Note that if a gene transitions n times, it is counted in n stage-pairs. Similarly, graph (H) shows the transitioning genes for the RNA-Seq dataset 
with c=0. 00025. The transitioning genes constitute 5% of all genes in that dataset. 45% of those genes transition in a single stage-pair. Of the 
remaining, 52% transition only in consecutive stage-pairs. 
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Figure '\3.Arabidopsis thaliana results using normalized expression levels. Graph (a) shows the hourglass score (normal and robust) as 
a function of the transition threshold c. Graph (b) shows the location of the hourglass waist (stage-pair) as a function of the transition threshold 
c. Graph (c) shows the Transcriptome Age Index of transitioning genes for five different values of c (chosen so that the number of genes with 
known age index assigned to each stage-pair is at least 290). Graph (d) shows the CDFs of the expression level absolute variations \S\ across 
successive stage-pairs. Note that in this case the hourglass waist (in terms of number of transitioning genes) appears in stage-pair (3,4), while 
the oldest genes appear in the next stage-pair. Graph (e) shows the transitioning genes with c=0.0001. The transitioning genes constitute 
7% of all genes in that dataset. 49% of those genes transition in a single stage-pair. Of the remaining, 36% transition only in consecutive 
stage-pairs. 



Figure 12a and Figure 12b show the hourglass resemblance score 
H (and its more robust variant) as function of c. Note that the H 
score is close to 1 for a wide range of c, confirming the presence of 
an hourglass-like structure in terms of the number of transitioning 
genes. Figure 12g and Figure 12h exhibit this pattern more clearly 
in the number of transitioning genes for a specific value of c. The 
two datasets also show reasonable agreement in terms of the assign- 
ment of transitioning genes to stage-pairs [Dataset. 10] (see Figure SI). 

Second, the location of the waist in this hourglass pattern, shown 
in Figure 12c and Figure 12d, occurs at the stage-pair (3,4) or 
(4,5), depending on c. This is roughly 8 hours after the formation 
of the zygote, and it includes the phylotypic stage for Drosophila 
melanogaster 9 . 

We have also estimated the evolutionary age of most of the tran- 
sitioning genes at each developmental stage-pair using the Tran- 
scriptome Age Index (TAI) metric 12 (see Methods). TAI is lower for 
older genes. Figure 12e and Figure 12f show the average TAI for 
transitioning genes, weighted by the expression level of each gene, 



at each stage-pair and for each dataset using three values of c. The 
TAI index follows the pattern that the model predicts, with older 
genes (lower TAI values) close to the waist of the hourglass. This 
result appears consistent with the main observation of Domazet- 
Loso and Tautz 12 , even though that study did not analyze transition- 
ing genes. 

The same approach was extended using transcriptome data 24 and 
TAI profiles 12 from Arabidops is. The results obtained from the anal- 
ysis (Figure 13) were consistent with the predictions of the model, 
as well as the results obtained from Drosophila data. 

Discussion 

Early studies of the developmental hourglass effect mostly ana- 
lyzed morphological and phenotypic similarities across species 2 27 . 
Recently, the focus has shifted towards genomic and molecular 
comparative studies 8 ' 9111219 that investigate conservation of gene 
expression variation, sequence conservation, selective constraint 
on coding sequences, and evolutionary gene "age". These stud- 
ies often report contradicting observations: some support strong 
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conservation in earlier developmental stages 5 619 , while others sup- 
port that strongest conservation occurs at a mid-developmental 
stage 715 . Nevertheless, the fact that the hourglass effect is observed 
in highly divergent species across deep phylogenetic scales (includ- 
ing fish, flies and plants), suggests that this observed pattern of con- 
servation may stem from fundamental organization principles. 

What these principles are has remained elusive. Earlier stages may 
be conserved because any changes therein could have large cascad- 
ing effects in later stages 5 28 29 . Later stages may experience less con- 
straint because as development progresses gene interactions become 
more modular, and so it is plausible that perturbations there have 
only local effects 1 . We refer to them as the "temporal" constraint 
model and the "spatial" constraint model, respectively, following 
Tian et al. 30 . 

In this paper, we developed an evolutionary model of development 
that combines some aspects of the previous two models. Regulatory 
perturbations at a certain stage can cause cascades of regulatory 
failures at subsequent stages (temporal model), while the likeli- 
hood that a gene regulates genes at a subsequent stage decreases 
as development progresses (spatial model). Our computational 
results lead to the following testable predictions: a) the number of 
transitioning genes during development follows an hourglass pat- 
tern, b) the evolutionary age of the transitioning genes also follows 
an hourglass pattern, with the oldest genes being at the waist of the 
hourglass, and c) the genes at the waist of that hourglass are the most 
essential, in the sense that their deletion maximizes the probability 
of developmental failure. We have relied on developmental gene 
expression profiles of Drosophila melanogaster and Arabidopsis 
thaliana to examine the predictions of the model. The analysis of 
that data agrees with the first two theoretical model predictions. The 
increased conservation of genes at the waist provides an indirect 
confirmation of the model's third prediction, regarding the essen- 
tiality of different genes. Further, our simulations confirm that the 
details of these regulatory perturbations, such as the probabilities of 
gene duplication and deletion or the parameter z in the regulatory 
failure probability, do not affect the results of the model, at least at 
the qualitative level. 

The use of DGENs in this work was only as an abstract tool to study 
the effect of gene regulatory perturbations in the developmental 
process. In future work, it is important to infer the actual DGEN of 
model organisms. This will require information about gene regula- 
tory interactions across time and space, but it should be possible for 
at least some developmentally well studied species 22 . Such DGENs 
would help to identify the specific genes that form the hourglass 
waist and their function. Additionally, an inferred DGEN would 
allow to directly test the increasing specificity assumption. 



Finally, we note that the hourglass effect (sometimes referred to as 
the "bow-tie" effect) has also been observed in other complex bio- 
logical and technological systems that exhibit hierarchical modu- 
larity and that are subject to evolutionary pressure or optimization 
tradeoffs 3134 . For instance, the Internet "protocol stack" is organ- 
ized in an hourglass structure 35 ; this pattern was not designed but 
it emerged through the competition between protocols that serve 
roughly the same function at each communication layer, during the 
last 30-40 years. In earlier work, we proposed an abstract model 
(EvoArch) that captures the evolution of protocol architectures and 
that predicts the emergence of an hourglass structure. Interestingly, 
both EvoArch and the model of this paper share the same princi- 
ple: the underlying hierarchical networks that control both systems 
should be increasingly sparser as complexity increases, i.e., the 
specificity of each complexity stage (or layer) should be increasing. 
In the future, we will further investigate this common organization 
principle between biological and technological systems. 

Data and software availability 

ZENODO: An evo-devo model for the developmental hourglass: 
code and data, doi: 10.528 l/zenodo.105 79 36 

License: GNU GPLv3 
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Figure S1. We evaluated the agreement between the two Drosophila datasets in terms of the transitioning genes assigned to each 
stage-pair, considering only those genes that appear in both datasets. Because the appropriate transition threshold may be different at 
each dataset, we use a different threshold for each dataset, say c 1 and c 2 . For each pair (c v c 2 ), we determine the transitioning genes at each 
stage-pair with the corresponding dataset (i.e., L - 1 pairs of gene sets), and then calculate the average Jaccard similarity across these L - 1 
pairs. The Jaccard similarity maps show this average across all stage-pairs for various threshold pairs (c v c 2 ). In graph (a) with normalized 
expression levels, when the two thresholds are roughly equal, the average Jaccard similarity is as high as 50%; this means that about 2/3 of 
the genes assigned to a certain stage-pair using one dataset are also assigned to the same stage-pair using the other dataset. Graph (b) 
shows a similar Jaccard similarity map for the case of absolute expression levels. [Dataset. 10]. 
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Figure S2. Drosophila results using absolute expression levels. Graphs (a) and (b) show the hourglass score (normal and robust) as a 
function of the transition threshold c for the two datasets. Graphs (c) and (d) show the location of the hourglass waist (stage-pair) as a function 
of the transition threshold c for the two datasets. Graphs (e) and (f) show the Transcriptome Age Index of transitioning genes for three different 
values of c (chosen so that the number of genes with known age index assigned to each stage is at least 10) for the two datasets. Graph (g) 
shows the transitioning genes for the Microarray dataset with c=5000. The transitioning genes constitute 21 % of all genes in that dataset. 62% 
of those genes transition in a single stage-pair. Of the remaining, 58% transition only in consecutive stage-pairs. Similarly, graph (h) shows the 
transitioning genes for the RNA-Seq dataset with c= 10000. The transitioning genes constitute 5% of all genes in that dataset. 45% of those 
genes transition in a single stage-pair. Of the remaining, 53% transition only in consecutive stage-pairs. [Dataset. 11]. 



Page 13 of 25 



FIOOOResearch 2014, 3:156 Last updated: 03 SEP 2014 



a) Normalized expression levels, Microarray 

■■, 



b) Normalized expression levels, RNA-Seq 




Stage 1 -2 
Stage 2-3 
Stage 3-4 
Stage 4-5 
Stage 5-6 
Stage 6-7 
Stage 7-8 
Stage 8-9 
Stage 9-10 



0.0001 



0.0002 



0.0003 



0.0004 



|S| 



c) Absolute expression levels, Microarray 



0.0005 




Stage 1 -2 
Stage 2-3 
Stage 3-4 
Stage 4-5 
Stage 5-6 
Stage 6-7 
Stage 7-8 
Stage 8-9 
Stage 9-10 



10000 



1 

0.9 
0.8 
0.7 
0.6 
0.5 
0.4 
0.3 
0.2 
0.1 
0 




Stage 1-2 
Stage 2-3 
Stage 3-4 
Stage 4-5 
Stage 5-6 
Stage 6-7 
Stage 7-8 
Stage 8-9 
Stage 9-10 



0.0001 



0.0002 0.0003 

|6| 



0.0004 



d) Absolute expression levels, RNA-Seq 



0.0005 




Stage 1-2 
Stage 2-3 
Stage 3-4 
Stage 4-5 
Stage 5-6 
Stage 6-7 
Stage 7-8 
Stage 8-9 
Stage 9-10 



10000 



|8| 



|8| 



Figure S3. Drosophila data: CDFs of the expression level absolute variations \3\ across successive stage-pairs, (a) normalized 
expressions, Microarray, (b) normalized expressions, RNA-Seq, (c) absolute expressions, Microarray, (d) absolute expressions, RNA-Seq. 
[Dataset. 12]. 
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Figure S4. Arabidopsis thaliana results using absolute expression levels. Graph (a) shows the hourglass score (normal and robust) as a 
function of the transition threshold c. Graph (b) shows the location of the hourglass waist (stage-pair) as a function of the transition threshold 
c. Graph (c) shows the Transcriptome Age Index of transitioning genes for five different values of c (chosen so that the number of genes with 
known age index assigned to each stage-pair is at least 1 70). Graph (d) shows the CDFs of the expression level absolute variations \S\ across 
successive stage-pairs. Graph (e) shows the transitioning genes with c=5000. The transitioning genes constitute 8% of all genes in that 
dataset. 48% of those genes transition in a single stage-pair. Of the remaining, 38% transition only in consecutive stage-pairs. [Dataset. 13]. 
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The manuscript by Akhshabi etal. proposes and describes in detail an evolutionary model that 
qualitatively reproduces the "Hourglass Effect"— a phenomenon (in the biological context) whereby there 
is increased morphological divergence at early and late stages of embryonic development separated by 
increased conservation in the phylotypic stage. 



The proposed model is based on a hierarchical network representation of the Gene Regulatory Network - 
the so-called Developmental Gene Execution Network (DGEN) - and evolution proceeds through a series 
of stochastic perturbations involving gene duplication, re-wiring and deletion. The manuscript suggests 
that the key factor that reproduces the hourglass effect is the monotonicity in "time" (increasing) of a 
function they call specificity s(i) (a measure of likelihood of a gene regulating descendant genes at later 
times) that competes against perturbative effects such that there is a "waist" of the number of transitions 
genes w(l) at intermediate times. They then test predictions of the model in two datasets of developmental 
data Drosophilia melanogaster and Arabidopsis Thaliana, finding good qualitative agreement. 

The manuscript is properly motivated and well written (with some caveats, see comments below). The 
model proposed is quite intuitive and compelling and represents an excellent first attempt in beginning to 
uncover the mechanism behind this intriguing phenomenon. I anticipate that this will have high impact 
(considering the hourglass effect transcends biological phenomena and is also found in "designed" 
systems such as the Internet protocol stack) in multiple fields. Consequently I strongly endorse this 
manuscript. 

I do however have some concerns about some of the details about the model and therefore propose the 
following that may improve the readability of the manuscript and the plausibility and robustness of the 
model: 

1 . The section describing the model is a bit dense, and can benefit from some rewriting. For example 
it is a bit hard to follow what the significance of a spatial domain is and its relevance to the model. 
In addition a bewildering plethora of parameters are introduced, that can be a bit overwhelming for 
the reader. Consequently I suggest putting all the important and relevant parameters in a table for 
easy access to the reader. The section on Rewiring (RW) is quite hard to read and will benefit from 
the introduction of explicit equations similar to that for P_{RF}. 

2. Which brings me to my second point about the assumptions of the model. It seems to me that two 
key factors lead to the hourglass effect: a) the montonocity of s(l) with stage I and likewise b) of 
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P_{RF} with r. While the choice of this dependence for the specificity is a simple linear 
dependence (which makes sense as a first pass), the choice of function for P_{RF} is highly 
non-trivial. Why not for example choose a simple linear dependence such as P_{RF} ~ r? The 
reason why I say this is that the authors would like to propose the most general and minimal model 
to explain this phenomenon. And it is not clear to me whether the specific choice of the functional 
dependences matter. For example with a non-linear choice of P_{RF} we see that a non-linear s(l) 
has the effect of shifting the waist. But what about a linear choice of P_{RF} and s(l)? Does that still 
preserve the hourglass effect? In my opinion a truly robust model will state that the main things that 
determine the hourglass are the monotonic dependencies in combination with non-linear or linear 
choices of the functions independent of their specific forms. This will allow the model to be applied 
to multiple settings. 

At the moment I'm a bit concerned that the non-trivial choice of P_{RF} 5 as it stands in the 
manuscript, might be key to the observed effect. If it is indeed so, then one must motivate why one 
must make this particular choice. Same goes for the slightly peculiar choice for the non-linearity of 
s(l) introduced later in the model. The authors will do well to explain in detail what motivated them 
to make such a choice. 

Additionally it would be helpful if the authors make it clearer (than they already have) that effects 
such as duplication and deletion are there only to make the model more realistic in terms of 
maintaining the number density of genes across different stages and play no part in the hourglass 
effect (or at least that's how I understand it). 

3. Some minor nitpicks: 

1 . In the section Model Description, the authors refer to the Wright-Fisher model. I was not 
aware of what this was and had to look it up, so a reference for the reader should be 
provided. 

2. The Hourglass score H is determined through Mann-Kendall statistics. A brief description of 
the method should be provided for the benefit of the reader (maybe in the supplementary 
material). Additionally, why this particular choice? Presumably any statistical test for 
monotonic trends should be robust to the parameters of the model. Would there be much of 
a difference is one used the Spearman rho for example? 

3. Finally I think the authors will do well to highlight the importance of Specificity right at the 
outset of the manuscript and move some simplified variant of the segment "What is the 
reason behind the hourglass shape of DGENS?" (page 7, line 1 1 , second column) to 
somewhere in the introduction. 

I have read this submission. I believe that I have an appropriate level of expertise to confirm that 
it is of an acceptable scientific standard, however I have significant reservations, as outlined 
above. 

Competing Interests: No competing interests were disclosed. 
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The article proposes and analyzes a network model that is capable of reproducing the hourglass effect in 
the number of genes that undergo a functional state transition as a function of the developmental stage. 
The hourglass effect refers to higher morphological divergence at earlier or later stages of embryonic 
developmental, compared to medium stages. The main idea of the model is to explain the waist of this 
hourglass (the number of "transitioning genes" w(l) having a minimum at medium developmental stages I) 
by an interplay between the impact of random rewiring events in the developmental gene execution 
network (DGEN) and stage specificity. The DGEN is a network representing regulatory interactions 
between genes at different developmental stages, while stage specificity s(l) defines the probability 1-s(l) 
with which a gene at stage I regulates a gene at stage 1+1 . One of the key assumptions in the model is that 
s(l) is an increasing function of I. As a result of the interplay between a decreasing impact of DGEN 
perturbations as a function of I, and increasing s(l), the lethality of perturbations is maximized at mid I's, 
leading to evolutionary incentives to minimize w(l) at those mid stages. The article consists of three parts. 
The first part describes the model in detail. The second part discusses extensive simulation results of the 
model. The third part presents an extensive analysis of existing developmental data on Drosophila 
melanogaster and Arabidopsis thaliana in the model context. 

The most interesting aspect of the article is that provides a possible intriguing explanation of the 
hourglass effect in developmental biology, which may foster future creative thinking and research in this 
important direction. The most obvious reservation one can express about the article is that the model 
cannot be currently refuted since it is formulated at the DGEN level, and there is currently no data from 
which a reliable DGEN reconstruction would be possible. The fact that some outcomes of model 
simulations qualitatively agree with available data does not directly validate either the model, or its main 
assumptions. The article does not contain such claims, however, that would not be supported by the data. 

My minor comments that may help to improve the article are: 

1 . The main idea of the model is buried somewhere on page 6 in the middle of a paragraph after 
sentence "What is the reason behind the hourglass shape of DGENs?" Why not explain this as 
early as possible, certainly in the Introduction, if not in the abstract? 

2. The third paragraph in the Introduction reads somewhat unclear and imprecise, as well as the 
"Developmental gene execution networks" section. 

3. When D(g) and Gamma notations first appear, it's not spelled out what they are. 

4. DGEN removals upon DF events have a clear meaning. But what do replacements correspond to 
in reality? 

5. The Mann-Kendall statistics could be briefly described, or at least a reference could be provided. 

I have read this submission. I believe that I have an appropriate level of expertise to confirm that 
it is of an acceptable scientific standard. 

Competing Interests: No competing interests were disclosed. 
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This study by Akhshabi etal. aims to provide the community with an explanatory theoretical model for the 
developmental hourglass phenomenon. Motivated by the question concerning the regulatory mechanism 
causing the developmental hourglass pattern during embryogenesis, the authors developed an 
evolutionary model of regulatory gene interactions during development to identify the conditions under 
which the molecular hourglass effect might emerge in general. 

For this purpose their evolutionary model focuses on hierarchical gene regulatory networks that control 
the corresponding developmental processes. Based on the output of the model the authors predict the 
emergence of an hourglass pattern in the structure of a temporal representation of the underlying gene 
regulatory network and correlate this effect with the evolution of protocol architecture found for the 
Internet's "protocol stack". 

Based on this universal finding the authors speculate that the developmental hourglass pattern might be a 
causal effect of a common organization principle during the establishment of complexity and is based on 
the underlying hierarchical network structure. 

The paper is very interesting, well written, and will potentially help the community to better 
understand the developmental hourglass phenomenon. However, we see some issues with the 
model assumptions and strongly suggest including the following aspects in a revised version. 

Summary of the points that have to be addressed: 

1 . Include a clear reference and paragraph to Raff's hypothesis and how it influenced the modeling of 
DGEN. 

2. Correlate and reference major assumptions of the DGEN model that seem plausible with findings 
from experimental studies that validate this plausibility. 

3. Correlate the output of the DGEN simulation with existing findings of studies addressing similar 
questions of the causality of the developmental hourglass phenomenon and discuss potential 
simplifications of the model assumptions that do not fit experimental results or that are not 
derivable via experiments yet. 

4. Include a clear statement how the exact values of the stage specificity (regulatory specificity) s(l) 
have been derived. 

5. State more clearly that s(l) is the crucial parameter determining the observed phenomena. 

6. Correct the TAI formula. 

7. Explain why merging OrthoDB results with OrthoMCL results might be plausible and why 
comparing it with results obtained from phylostratigraphy is reliable. 

8. Discuss more clearly how the DGEN model can help to find causal processes of the 
developmental hourglass phenomenon based on recent knowledge: the three predictions of the 
DGEN model a) - c) mentioned in the Discussion section seem not integrative enough to allow 
future studies to correlate existing results provided by the community with the predictions obtained 
by the DGEN model. 
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Specific comments to the separate sections: 

In general, most of our criticisms concern the model description and the methods section. Both are 
integral to the validity of the obtained results, which are nicely reported and discussed. 

Introduction 

Modeling the evolution of embryonic development as directed acyclic graph (DAG) in which the nodes 
correspond to state-transitioning genes and the edges model directed regulatory interactions causing 
significant activity change at the corresponding stage seems plausible. The authors refer to this model as 
DGEN. A testable hypothesis proposed by the authors is that if regulatory genes become increasingly 
function-specific as development progresses, the network gradually should become sparser in that 
direction and evolutionary older genes should be concentrated at the waist of the hourglass (Figure 1). 

This hypothesis is essential to detect possible regulatory mechanisms that cause the developmental 
hourglass pattern. We therefore suggest including references to existing biological studies investigating 
the phenomenon of morphological complexity being correlated with regulatory specialization to support 
the validity and plausibility of their DGEN model assumptions. 

A biological motivation of this hypothesis comes from Rudolf A. Raff who proposed that early 
embryogenesis is governed by global pattern formation processes and that during later development, the 
embryo is being organized into developmental modules, allowing each module to differentiate 
autonomously. Raff furthermore proposed that the phylotypic stage marks the transition from early global 
development to the modular mode later in development 1 . 

The authors do cite Raff in the first sentence of their introduction, but should state more clearly in which 
manner their model was inspired by Raff's hypothesis. 

Numerous studies investigated the phenomenon that increasing morphological complexity is correlated 
with increasing modularity and specialization of the underlying regulatory apparatus 2 . Here the authors 
predict that the main condition leading to the appearance of the hourglass pattern is that the DGEN 
should gradually get sparser as development progresses, with general-purpose regulatory genes at the 
earlier developmental stages and highly specialized regulatory genes at the later stages, expressing the 
oldest genes during mid development. 

We suggest to also include findings investigating the functionality and age of the genes that have been 
shown to be active during mid embryogenesis to correlate potential outcomes of the DGEN prediction that 
evolutionary older genes should be concentrated at the waist of the hourglass with results obtained by 
recent studies 34567891011 12 jh ese aforementioned studies were able to show that during mid 
embryogenesis evolutionary old genes or evolutionary old processes that are shared within and between 
phyla are most active. 

In our opinion it would be very interesting to test the degree of intersection between genes predicted to be 
causal for the waist of the DGEN hourglass and evolutionary conserved genes that have been found to be 
expressed during mid embryogenesis. We are aware that the model uses the number of genes at each 
developmental stage combined with the age of the corresponding genes to observe constraints during 
development. Hence, only the authors can estimate how much effort would be needed to identify such 
intersecting genes. 
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If the search for intersecting genes returned by the DGEN model during mid development and genes 
found by previous studies should be possible, the authors could integrate a motivation section to find 
these intersecting genes within the Introduction section. 

Model description 

The model description is quite intuitive. However, as the validity of the whole study depends on the model 
parameters chosen, several points need clarification including the rationale and motivation to choose 
exactly these parameters. 

One of the most crucial parameters of the model is the regulatory specificity s(l), with 0 <= s(l) <= 1 . We 
think, that adding a more detailed paragraph discussing how s(l) is being estimated, would be highly 
beneficial due to the importance for the parameter s(l) for the entire model and predictions. In Model-1 -4 
the authors choose three different values for s(l): In Model-1 : s(l) = 0.5 and in Model-2-4: s(l) = l/L, and in 
Figure 9 a non-linear function s(l) = 0.9 - (0.8 / 1 + exp(\gamma - 1) ). A gene g at stage I is modeled to act 
as upstream regulator for a gene g' at stage I + 1 with probability s'(l) = 1 - s(l). From this follows that each 
potential upstream regulator in stage I has the same probability s'(l) to regulate genes g' in stage I + 1 . In 
other words, Model-1 assigns each potential upstream regulator with the same probability s'(l) = 1 - 0.5 = 
0.5 to regulate genes g' in stage I + 1 , in Model-2: s'(l) = 1 - l/L = const., and for the non-linear function s(l) 
: s'(l) = 1 - (0.9 - (0.8 / 1 + exp(\gamma - 1) )) = const. 

It would be advantageous, if the authors could include motivations for the three cut-offs, especially how 
they derived the sigmoid-like mathematical function (Figure 10, legend caption) they later denote as 
non-linear function s(l). We do see that modeling the stage specificity with a constant value, an increasing 
value, and a non-linear function is plausible, but a clear motivation or derivation of the exact values would 
be beneficial for a better understanding of the modeling process. 

For the probabilities P DL and P DP the explicit formula to compute the corresponding probabilities could 
also be included within the paragraph analogous to P RF . This would enable a clearer reproducibility. 

The authors discuss that their major assumption is that the regulatory specificity increases substantially 
as development progresses, hence the DGEN becomes gradually sparser along the developmental time 
axis starting with s(1) ~ 0, ... , s(L) ~ 1 . Here the authors should include references to experimental 
observations that support their assumption. Since Models 2-4 rely on this assumption, published 
experimental studies in line with this assumption need to be referenced. Otherwise, the models applied 
would seem largely academic. 

Methods 

Hourglass score H 

Here w(l) is defined as the number of transitioning genes in stage I. In Figure 3 and in section 'Hourglass 
shape' w(l) is referred to as stage width. In case the definition of w(i) is equivalent to the stage width, it 
would be beneficial for the reader to see the term stage width in the same sentence as the initial definition 
of w(l). A correct understanding of the stage width is important to understand the measurement of H. 

Furthermore, the univariate Mann-Kendall statistic should be referenced to allow a non-statistical 
community to better understand this way of statistical modeling. 
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In this section the authors describe and perform two different methods to assign each gene a 
corresponding evolutionary age. For D. melanogaster Vney collected groups of orthologs from OrthoDB 
and OrthoMCL, whereas for A thaliana they obtained phylostratum assignments from Quint et al. 10 that 
was derived by phylostratigraphy[ref-13.]We do not understand the motivation of merging orthologous 
gene groups obtained from OrthoDB and OrthoMCL, as both approaches apply different parameters and 
heuristics to determine orthologous genes. Was this done to maximize the number of potential orthologs 
to be included in the study? In addition, why did the authors use a different approach for the plant gene 
set? This seems somewhat inconsistent and statement concerning the plausibility to compare results 
returned by the DGEN model based on age assignments obtained from phylostratigraphy (A thaliana) 
with OrthoDB and OrthoMCL merged orthologs (D. melanogaster) would be appreciated. 

Furthermore, the database specific parameters they used in OrthoDB and OrthoMCL should be 
mentioned as well as the database version they used, e.g. OrthoDB2, or OrthoDB3, ... , or OrthoDB6. 

Age index for each stage-pair 

There seems to be a typo within the TAI formula. As defined by Domazet-Loso and Tautz 6 and adapted to 
the author's notation the formula of the TAI should begin TAI(I) not TAI(1). 

Results 

The authors show based on their findings returned by Model-2-4 that the increasing specificity 
assumption is the key property behind the developmental hourglass effect. We think that this statement is 
not clear enough. The authors could phrase more clearly what they mean by "key property behind the 
hourglass effect" inferred from Model-2-4. 

Age of genes 

Here the authors motivate that defining the age of a gene as A(g) = i-t 0 (g) is plausible, because a 
rewiring event may give a gene a new function, at least in terms of its upstream and downstream 
regulators. This statement needs validation from the literature and from experimental studies that already 
tested this hypothesis for specific rewiring events. 

Furthermore, the authors show (Figure 8c) that the evolutionary age at stage I follows the same pattern as 
the lethality probability: gradually increasing until the waist of the hourglass, and then gradually 
decreasing. This finding should be correlated with the findings by Galis and Metz 3 and should be 
discussed more clearly in the Discussion. 

In the first sentence of the last paragraph before the Discussion, "... and TAI profiles from Arabidopsis" 
should be referenced to reference 1 1 and not 12. 

Discussion 

The authors discuss that the observed pattern of conservation (developmental hourglass) may stem from 
fundamental organization principles and that the exact origin of these principles remains elusive. Here the 
authors should include existing studies that also aim at predicting a universal organization principle during 
development (modularization, etc.) and what exact problems are preventing studies to elucidate this 
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universal phenomenon. We do note, that the authors cite studies that aim to explain the effects causing 
the early phase of the developmental hourglass 141516 and studies aiming to explain effects causing the 
late phase of the developmental hourglass 1 , but we think that it should be explicitly stated that the model 
is designed to test aspects of Raffs initial hypothesis and that also more recent studies have been 
investigating the underlying reasons for the pattern in the early phase of the hourglass 17 . 
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We have read this submission. We believe that we have an appropriate level of expertise to 
confirm that it is of an acceptable scientific standard, however we have significant reservations, 
as outlined above. 
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